##Replication File for "Authoritarian Audiences, Rhetoric, and Propaganda in International Crises: Evidence from China"
##Authors: Jessica Chen Weiss and Allan Dafoe

##Replication Real History
rm(list=ls())

#Packages
library(ggplot2)
require(reshape2)
require(coefplot)
library(robustbase)

#hyp.Rdata/real.Rdata cannot be shared due to privacy concerns, instead the de-identified datasets with simulated demographics (hyp_toshare and real_toshare) can be used to replicate the paper's results.

setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")



#------------------------------------------------------------------------------------------
## Constructing Publishable dataset. 
#------------------------------------------------------------------------------------------
#load("./MASTER CODE/real.Rdata")

#keepvars <- c("asc", "asc0", "pre.questions", "asc.or", "his", "pro", "ADIZ", "ADIZp", "eli.f",
#            "eli.c", "asc0.v2", "na1.v2", "na2.v2", "na3.v2", "na2.v.dn", "na3.v.dn", "start.time",
#             "start.time.n", "start.time.n2", "start.time.n3", "start.time.swd",
#              "na1.v", "na2.v", "na3.v", "att", "rac",  "rac.or", "ctri", 
#              "V1") ##variables that are not identifiable

#data.share <- d[, names(d) %in% keepvars] #keeping only shareable data

#write.csv(data.share, "./MASTER CODE/real_toshare.csv")

d <- read.csv("./MASTER CODE/real_toshare.csv")

#analysis data keeps only non-identifiable variables. 
#generating random numbers of identifiable variables
set.seed(2018)

d$gender <- sample(c(0,1), size=nrow(d), replace=TRUE, prob=c(.37,.63))
d$educ <- sample(seq(from = 1, to = 7, by = 1), size = nrow(d), replace = TRUE, prob=c(.0013, .0007, .417, .071, .452, .050, .007))
d$age <- sample(c(11,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33 , 34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69, 70,72,73,76,80,82,85),
                size=nrow(d), replace=TRUE, prob=c(0.0002, 0.0005, 0.0052, 0.0087, 0.0118, 0.0124, 0.0153, 0.0148, 0.0156, 0.0319, 0.0207, 0.0266, 0.0284, 0.0240, 
                                                   0.4316, 0.0237, 0.0242, 0.0286, 0.0205, 0.0328, 0.0182, 0.0183, 0.0148, 0.0158, 0.0163, 0.0136, 0.0126, 0.0126,
                                                   0.0111, 0.0139, 0.0082, 0.0067, 0.0064, 0.0062, 0.0040, 0.0045, 0.0050, 0.0029, 0.0015, 0.0020, 0.0032, 0.0030,
                                                   0.0024, 0.0025, 0.0020, 0.0018, 0.0012, 0.0027, 0.0024, 0.0015, 0.0015, 0.0007, 0.0007, 0.0003, 0.0007, 0.0002,
                                                   0.0003, 0.0002, 0.0002, 0.0002, 0.0002))
d$partner <- sample(c("A", "B"), size=nrow(d), replace=TRUE, prob=c(.295, .705)) 


#Gen other variables 
d$gender.o <- d$gender
d$gender.m <- is.na(d$gender)
d$gender[d$gender.m==TRUE] <- 1

d$age.o <- d$age
d$age.m <- is.na(d$age)
d$age[d$age.m==TRUE] <- 30

d$educ.o <- d$educ
d$educ.m <- is.na(d$educ)
d$educ[d$educ.m==TRUE] <- 3


write.csv(d, "./MASTER CODE/realanalysis.csv")
save(d, file="./MASTER CODE/realanalysis.Rdata")


#-----------------------------------------------------------------------------------------
#Replication Main Text and Appendix
#-----------------------------------------------------------------------------------------
{{format(Sys.Date(), format="%B %d %Y")}}

#+ opts, include=FALSE,eval=TRUE
knitr::opts_chunk$set(eval=TRUE, echo=TRUE, error=TRUE,message=FALSE, warning=FALSE,
                      fig.width=15, fig.height=5,width=60,dev="pdf"
)

rm(list = ls())

setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

d<-read.csv("./MASTER CODE/realanalysis.csv")
  rb <- rep(1, length(d[,1]))

#Figures location
fig.loc <- "../../Figures/"

#Packages
require(ggplot2)
require(reshape2)
require(coefplot)
require(robustbase)
require(rms)
require(sandwich)
library(lubridate)
library(stargazer)

w <- 5
h <- 4

#finding value of Z that gives 10% one sided test
xs <- seq(-2,0, 0.01)
xs[which(pnorm(xs)-0.05==min(abs(pnorm(xs)-0.05)))]
pnorm(-1.64)

innerCI.n <- 1.64 #set so exclusion of CI implies two sided significance of 0.1
outerCI.n <- 1.96 #set so exclusion of CI implies two sided significance of 0.05



# Customized Stargazer Function, without asterisk inflation
stargazerAD <- function(table, title="Default", filename="Figures/default.tex", dep.var.labels=NULL) {
  a <- stargazer(table, title=title,  out=filename, 
                 align=TRUE, star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c("\\dagger", "*", "**", "***"), 
                 notes="$^{\\dagger} p<0.1$; $^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$", notes.append=FALSE, single.row=TRUE)
  return(a)
}

results.fun <- function(data, r, model, title.w, fig.loc, coefficients=NULL, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=10, h=6) {
  d <- data
  rob.results <- robcov(ols(model, dat=d[r,], x=TRUE))
  m1 <-  lm(model,  dat=d[r,])
  homo.results <- summary(m1)
  n <- length(d[r,1])
  title <- paste(title.w, ", n=", n, sep="")
  pf <- coefplot(m1, title=title, 
                 intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n,
                 coefficients=coefficients) + theme_bw()
  fig.results <- pf
  list.results <- list("rob.results"=rob.results, "homo.results"=homo.results, "title"=title, "fig.results"=fig.results, "model"=m1, m1)
  return(list.results)
}



## Number of observations ##
n.obs <- sum(table(d$asc))
cat(n.obs, file=paste(fig.loc,"Hist-n.tex", sep=""))

n.obs.att <- sum(d$att) 
cat(n.obs.att, file=paste(fig.loc,"Hist-n-att.tex", sep=""))



#From PAP: 
#For our primary specification we will estimate the effects of our manipulations using OLS with heteroscedasticity consis- tent standard errors, conditioning on the other experimental conditions to increase power. Hypothesis tests will be one-sided when appropriate. We will also report our estimates of the average causal effect when not conditioning on anything else, an estimand termed the average marginal component effect by Hainmueller, Hopkins, and Yamamoto.37
#Our primary test of these hypotheses will involve OLS regression, controlling for the other experimental conditions, and in the real-history design for the pre-scenario answers to the same questions. 
#robcov(ols(y ~ a.1 * a.2, toyData, x = TRUE)) ## heteroskedastic (ATE)

### Results for Approval #### 
#' ### Results for Approval #### 


#r <-  rep(1, length(d[,1]))
r <- !is.na(d$asc) & rb #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Hist"
outcome <- "Approval"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")
coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- asc ~ pre.questions + asc.or + his + pro + ADIZ + ADIZp + eli.f + eli.c
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()
m11 <- results[5]


pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()

pro.results <- coef(summary(results$model))["pro",]
pro.results <- rbind(pro.results, coef(summary(results$model))["ADIZp",])
pro.results <- as.data.frame(pro.results)
names(pro.results) <- c("est", "se", "t", "p")
pro.results$var <- c("EP-3", "ADIZp")
pro.results$dv <- "approval"
pro.results$cov <- "no covariates"

toprint <- paste("p_{+}=",round(pro.results[1,4]/2,3),  sep="")
write(toprint, file=paste(fig.loc,type,outcome, "EP3", "p",".tex", sep=""))

statement.results <- coef(summary(results$model))[c("ADIZ", "his", "eli.f", "eli.c"),]


#Note: effect of ADIZ and ADIZp become positive after Nov 2015, consistent with any tough action being approved of.
#Note, hypothesis tests should be read off of robcov function, to allow for heteroscedasticity. For convenience, the coefficient plot shows plots confidence intervals based on assuming homoscedasticity. 




### Now controlling for pre.questions, partner, time, and demographics ####
#' ### Now controlling for pre.questions, partner, time, and demographics ###
r <- !is.na(d$asc) & rb  #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Hist"
outcome <- "Approval"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")
coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- asc ~ his + pro + ADIZ + ADIZp + eli.f + eli.c +
  partner + pre.questions + asc.or + asc0.v2  + na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()
m12 <- results[5]

stargazerAD(c(m11, m12), title="Effect on Approval, History", filename=paste(fig.loc,title.f,"-table",".tex", sep=""))



smalld <- d[,c("asc", "asc0", "his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c", 
               "partner", "asc.or", 
               "na1.v", "na2.v", "na3.v", "na2.v.dn", "na3.v.dn", "gender.o", "educ.o", "age.o",
               "age.m", "gender.m", "educ.m")]
#, "start.time.n", "start.time.n2", "start.time.n3", "start.time.swd"
stargazer(smalld, title="Summary Statistics", out=paste(fig.loc, "summary-statistics-history",".tex", sep=""))



pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()

pro.results <- rbind(pro.results, 
                     c(as.numeric(coef(summary(results$model))["pro",]), "EP-3", "approval", "covariates"))
pro.results <- rbind(pro.results, 
                     c(as.numeric(coef(summary(results$model))["ADIZp",]), "ADIZp", "approval", "covariates"))

statement.results <- rbind(statement.results, coef(summary(results$model))[c("ADIZ", "his", "eli.f", "eli.c"),])
row.names(statement.results) <- paste(row.names(statement.results), rep(", ", 4), c(rep("approval",4), rep("approval, cov", 4)), sep="")

##Analysis replicated among only those who passed attention check
#r <-  rep(1, length(d[,1]))
r <- !is.na(d$asc) & rb #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Inattentive-Hist"
outcome <- "Approval"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")
coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- asc ~ pre.questions + asc.or + his + pro + ADIZ + ADIZp + eli.f + eli.c
results <- results.fun(data=d[d$att==1,], r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
m11_inattentive <- results[5]



##Analysis replicated among only those who passed attention check
### Now controlling for pre.questions, partner, time, and demographics ####
#' ### Now controlling for pre.questions, partner, time, and demographics ###
r <- !is.na(d$asc) & rb  #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Inattentive-Hist"
outcome <- "Approval"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")
coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- asc ~ his + pro + ADIZ + ADIZp + eli.f + eli.c +
  partner + pre.questions + asc.or + asc0.v2  + na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd
results <- results.fun(data=d[d$att==1,], r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
m12_inattentive <- results[5]
stargazerAD(c(m11_inattentive, m12_inattentive), title="Effect on Approval, History (Inattentive respondents omitted)", filename=paste(fig.loc,title.f,"-table",".tex", sep=""))



#analysis of distribution of inattentiveness 


r <- !is.na(d$asc) & rb #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Random-Inattentive-Hist"
outcome <- "Approval"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")
coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- att ~ pre.questions + asc.or + his + pro + ADIZ + ADIZp + eli.f + eli.c
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
m31_inattentive <- results[5]



### Now controlling for pre.questions, partner, time, and demographics ####
#' ### Now controlling for pre.questions, partner, time, and demographics ###
r <- !is.na(d$asc) & rb  #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Random-Inattentive-Hist"
outcome <- "Approval"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")
coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- att ~ his + pro + ADIZ + ADIZp + eli.f + eli.c +
  partner + pre.questions + asc.or + asc0.v2  + na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn + gender + educ + age +
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2] 
results$fig.results
m32_inattentive <- results[5]
stargazerAD(c(m31_inattentive), title="Effect on Attentiveness, History", filename=paste(fig.loc,title.f,"-table",".tex", sep=""))

### Analyzing Resolve ####
#' ### Analyzing Resolve ###
r <- !is.na(d$rac) & rb #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Hist"
outcome <- "Resolve"
cont <- ""
title.w <- paste(type, ", Effect on ", outcome, cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")
coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- rac ~ pre.questions + rac.or + his + pro + ADIZ + ADIZp + eli.f + eli.c
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2]
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()

pro.results <- rbind(pro.results, 
                     c(as.numeric(coef(summary(results$model))["pro",]), "EP-3", "resolve", "no covariates"))
pro.results <- rbind(pro.results, 
                     c(as.numeric(coef(summary(results$model))["ADIZp",]), "ADIZp", "resolve", "no covariates"))

toprint <- paste("p_{+}=",round(as.numeric(pro.results[5,4])/2,3),  sep="")
write(toprint, file=paste(fig.loc,type,outcome, "EP3", "p",".tex", sep=""))

toprint <- paste("p_{+}=",round(as.numeric(pro.results[6,4])/2,3),  sep="")
write(toprint, file=paste(fig.loc,type,outcome, "ADIZp", "p",".tex", sep=""))


statement.results <- rbind(statement.results, coef(summary(results$model))[c("ADIZ", "his", "eli.f", "eli.c"),])
r <- 9:12
row.names(statement.results)[r] <- paste(row.names(statement.results)[r], rep(", ", 4), c(rep("resolve",4)), sep="")



r <- !is.na(d$rac) & rb #& d$start.time > "2015-11-04 21:02:13 EST"
type <- "Hist"
outcome <- "Resolve"
cont <- "Covariates"
title.w <- paste(type, ", Effect on ", outcome, ", Controlling ", cont, sep="")
title.f <- paste(type, "-", outcome, cont, sep="")

coefficients.list<- c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c")
model <- rac ~ pre.questions + his + pro + ADIZ + ADIZp + eli.f + eli.c +
  partner + pre.questions + rac0.v2 + rac.or + na1.v2 + na2.v2 + na3.v2 + na2.v.dn + na3.v.dn +
  gender + educ + age + age.m + gender.m + educ.m + 
  start.time.n + start.time.n2 + start.time.n3 + start.time.swd
results <- results.fun(data=d, r=r, model=model, title.w=title.w, fig.loc=fig.loc, coefficients=coefficients.list, innerCI.n=innerCI.n, outerCI.n=outerCI.n, w=w, h=h)
results[1]
results[2]
results$fig.results
coefplot(lm(model, data=d[r,]), title=title.w, 
         intercept=FALSE, lwdOuter=0.5, innerCI=innerCI.n, outerCI= outerCI.n, data=d) + theme_bw()

pdf(paste(fig.loc,title.f,".pdf", sep=""), width=w, height=h)
results$fig.results
dev.off()

pro.results <- rbind(pro.results, 
                     c(as.numeric(coef(summary(results$model))["pro",]), "EP-3", "resolve", "covariates"))
pro.results <- rbind(pro.results, 
                     c(as.numeric(coef(summary(results$model))["ADIZp",]), "ADIZp", "resolve", "covariates"))


toprint <- paste("p_{+}=",round(as.numeric(pro.results[7,4])/2,3),  sep="")
write(toprint, file=paste(fig.loc,type,outcome, "EP3", "p", "Cov", ".tex", sep=""))

toprint <- paste("p_{+}=",round(as.numeric(pro.results[8,4])/2,3),  sep="")
write(toprint, file=paste(fig.loc,type,outcome, "ADIZp", "p","Cov",".tex", sep=""))

statement.results <- rbind(statement.results, coef(summary(results$model))[c("ADIZ", "his", "eli.f", "eli.c"),])
r <- 13:16
row.names(statement.results)[r] <- paste(row.names(statement.results)[r], rep(", ", 4), c(rep("resolve, cov",4)), sep="")


## Distribution of Approval ##

hist(d$asc0, main="Historical Approval (Pre)", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE, col="blue", xaxt='n', xlab="", cex=0.8)
axis(side=1, at=c(1, 3, 5), labels=c("Strongly \n disapprove", "Neither approve \n nor disapprove", "Strongly \n approve"))

pdf(paste(fig.loc,"Hist-Approval(Pre)-Histogram",".pdf", sep=""), width=w, height=h*1.6)
par(mar=c(3.1,4.1,4,2.1))
hist(d$asc0, main="Approval (Pre) in History", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE, col="blue", xaxt='n', xlab="", cex=0.8)
axis(side=1, at=c(1, 3, 5), labels=c("Strongly \n disapprove", "Neither approve \n nor disapprove", "Strongly \n approve"))
dev.off()



hist(d$asc, main="Historical Approval (Post)", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE,  col="blue", xaxt='n', xlab="")
axis(side=1, at=c(1, 3, 5), labels=c("Strongly disapprove", "Neither approve nor disapprove", "Strongly approve"))

pdf(paste(fig.loc,"Hist-Approval(Post)-Histogram",".pdf", sep=""), width=w, height=h*1.6)
par(mar=c(3.1,4.1,3,2.1))
hist(d$asc, main="Approval (Post) in History", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE, col="blue", xaxt='n', xlab="", cex=0.6)
axis(side=1, at=c(1, 3, 5), labels=c("Strongly \n disapprove", "Neither approve \n nor disapprove", "Strongly \n approve"))
dev.off()



# #par(oma=c(0,0,-2,0))
# title("Historical Approval (Post)", line = 2)

#How many respondents got ADIZp
sum(d$ADIZp==1 & d$att==1)
sum(d$ADIZ==1 & d$att==1)
table(d$ADIZ, d$ADIZp)
sum(d$pro==1 & d$att==1)


pro.results$var <- as.factor(pro.results$var)
pro.results$dv <- as.factor(pro.results$dv)
pro.results$cov <- as.factor(pro.results$cov)
pro.results$est <- as.numeric(pro.results$est)
pro.results$se <- as.numeric(pro.results$se)
pro.results$t <- as.numeric(pro.results$t)
pro.results$p <- as.numeric(pro.results$p)
pro.results$names2 <- paste(pro.results$var," \n ", pro.results$cov, sep="")
pro.results


statement.results <- as.data.frame(statement.results)
statement.results$names <- as.factor(row.names(statement.results))
statement.results$names2 <- as.factor(rep(c(rep("No Covariates",4),rep("Covariates",4)),2))
statement.results$names3 <- as.factor(rep(c("Vague Threat", "Nationalist History", "Biding Time", "Cost of War"),2))
statement.results$names3 <- factor(statement.results$names3, 
                                   levels=rev(c("Vague Threat", "Nationalist History", "Cost of War", "Biding Time")))
head(statement.results)
# factor(x$name, levels = x$name[order(x$val)])

statement.results$dv <- c(rep(0,8), rep(1,8))
names(statement.results) <- c("est", "se", "t", "p", "names", "names2","names3", "dv")
statement.results$names4 <- paste(statement.results$names3, "\n ", statement.results$names2, sep="")
statement.results


pvalues <- -1

r <- !is.na(d$asc) & rb 
n <- length(d[r,1])

#statements, statements cov, mob, mob covs, adiz, adiz covs, 

type <- "Hist"
r2 <- statement.results$dv==0 & statement.results$names3=="Vague Threat"   #& statement.results$names2=="No Covariates" 
dv.n <- "Approval"
title.w <- "Effect on"
covs <- ""
main.t <- "1StateA-"
ymin.n <- -1; ymax.n <- .4

source("./MASTER CODE/6-coef-figuresb.R", echo=TRUE)


type <- "Hist"
r2 <- statement.results$dv==0 & statement.results$names3!="Vague Threat"   #& statement.results$names2=="No Covariates" 
dv.n <- "Approval"
title.w <- "Effect on"
covs <- ""
main.t <- "1StateB-"
ymin.n <- -1; ymax.n <- .4
source("./MASTER CODE/6-coef-figuresb.R", echo=TRUE)



type <- "Hist"
r2 <- statement.results$dv==0 & statement.results$names2=="Covariates"
dv.n <- "Approval"
title.w <- "Effect on"
covs <- "Covs"
main.t <- "1State-"
ymin.n <- -1; ymax.n <- .4
source("./MASTER CODE/6-coef-figuresb.R", echo=TRUE)


#Demographics 


rm(list = ls())

setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

d <- read.csv("./MASTER CODE/realanalysis.csv")

#Figures location
fig.loc <- "../../Figures/"

w <- 5
h <- 4
innerCI.n <- 1.64 #set so exclusion of CI implies two sided significance of 0.1
outerCI.n <- 1.96 #set so exclusion of CI implies two sided significance of 0.05



shadowtext <- function(x, y=NULL, labels, col='white', bg='black',
                       theta= seq(pi/4, 2*pi, length.out=8), r=0.1, ... ) {
  
  xy <- xy.coords(x,y)
  xo <- r*strwidth('A')
  yo <- r*strheight('A')
  for (i in theta) {
    text( xy$x + cos(i)*xo, xy$y + sin(i)*yo, 
          labels, col=bg, ... )
  }
  text(xy$x, xy$y, labels, col=col, ... )
}


stat_sum_df <- function(fun, geom="crossbar", ...) {
  stat_summary(fun.data=fun, colour="red", geom=geom, width=0.2, ...)
}


pdf(paste(fig.loc,"Hist-Approval(Pre)-Hist",".pdf", sep=""), width=w, height=h)
par(mar=c(3.1,4.1,4,2.1))
hist(d$asc0, main="Approval (Pre) in Hist", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE, col="blue", xaxt='n', xlab="", ylim=c(0,0.5))
axis(side=1, at=c(1, 3, 5), labels=c("Strongly \n disapprove", "Neither approve \n nor disapprove", "Strongly \n approve"))
dev.off()


plot(density(d$rac0[d$pre.questions==1]), main="Resolve")


#Second, we would expect to see high levels of agreement with the question $cc$: To what extent do you agree or disagree with the following statement: ``Criticism of government policy is unhelpful.'' To the extent that we do not see see overwhelming high levels of $cc$, we gain confidence that respondents are not self-censoring. 
hist(as.numeric(as.character(d$cc)))
hist(as.numeric(as.character(d$ccr)))
summary(lm(d$ccc ~ d$cc.or)) #A small proportion of the answer can be attributed to the order of the answer options. 

pdf(paste(fig.loc,"Hist-Criticism",".pdf", sep=""), width=w, height=h)
par(mar=c(3.1,4.1,4,2.1))
hist(d$ccc, main="Criticism of the Government is Unhelpful", breaks=seq(0.5, 5.5, 1), freq=FALSE, labels=TRUE, col="blue", xaxt='n', xlab="", ylim=c(0,0.5))
axis(side=1, at=c(1, 3, 5), labels=c("Strongly \n disagree", "Neither disagree \n nor agree", "Strongly \n agree"))
dev.off()


#Pre Approval doesn't correlate with support of government criticism
#summary(lm(d$ccc ~ d$asc0)) 
#Concern for national honor slightly negatively associated with support of government criticism
#summary(lm(d$ccc ~ d$na1.v)) 

## Number of observations who passed attention checks
sum(d$att)
sum((d$ctri <= 2) %in% TRUE)

hist(d$age)
age <- d$age

#65% male of attentives; 75% male of initial
prop.table(table(d$gender))

d$educ.f <- as.factor(d$educ)
d$educ.f[d$educ.m==1] <- NA
library(plyr)
# 01 No formal education 02 Elementary school 03 Middle school
# 04 High school
# 05 College 06 Masters 07 Doctoral
d$educ.f <- revalue(d$educ.f, c("1"="No Educ", "2"="Elementary", "3"="Middle", "4"="High School", "5"="College", "6"="MA", "7"="Doctoral"))

prop.q.data <- table(d$educ.f)/sum(!is.na(d$educ.f))
prop.ci.data <- c(0,.119,.36,.312,.209, 0, 0)
prop.data <- rbind(prop.q.data, prop.ci.data)
row.names(prop.data) <- c("Our Sample", "Chinese \nInternet Users")
pdf(paste(fig.loc,"hist_education.pdf", sep=""), width=8, height=6)
barplot(prop.q.data, main="Education", cex.names=0.7)
dev.off()



pdf(paste(fig.loc,"hist_education_comp.pdf", sep=""), width=8, height=6)
barplot(prop.data, main="Education", cex.names=0.7, beside=TRUE, legend=rownames(prop.data))
dev.off()
p <- recordPlot()
p



d$educ.f2 <- d$educ.f
d$educ.f2[d$educ.f=="MA" | d$educ.f=="Doctoral"] <- "MA"
d$educ.f2[d$educ.f=="No Educ"] <- "Elementary"
d$educ.f2 <- revalue(d$educ.f2, c("MA"="Graduate"))
d$educ.f2 <- droplevels(d$educ.f2, c("No Educ","Doctoral"))
prop.q.data <- table(d$educ.f2)/sum(!is.na(d$educ.f2))

#Primary, Junior, Senior, 2 year college + 4 year college, Graduate
prop.huang.data <- c(0.006, 0.038, 0.103, 0.307 + 0.507, 0.04)
prop.data <- rbind(prop.q.data, prop.huang.data)
row.names(prop.data) <- c("Our Sample", "Huang Data")

barplot(prop.data, main="Education", cex.names=0.7, beside=TRUE, legend.text=rownames(prop.data), args.legend=list(x=6, y=0.8))

pdf(paste(fig.loc,"hist_education_comp_Huang.pdf", sep=""), width=8, height=6)
barplot(prop.data, main="Education", cex.names=0.7, beside=TRUE, legend.text=rownames(prop.data), args.legend=list(x=6, y=0.8))
dev.off()


#### Age analysis ####
r <- d$age.m==0
pdf(paste(fig.loc,"hist_age_den.pdf", sep=""), width=8, height=6)
den <- density(age[r]) # returns the density data 
plot(den, col="cyan3", main="Age of Respondents") 
polygon(den, col="cyan3", border="black")
rug(jitter(age, amount=0.5))
dev.off()

##Actual Chinese age distribution based on 2010 census, eye-ball estimated from this figure:
#http://en.wikipedia.org/wiki/File:China_Sex_By_Age_2010_census.png

#From: https://www.cia.gov/library/publications/the-world-factbook/fields/2010.html
#0-14 years: 17.1% (male 124,340,516/female 107,287,324) 
#15-24 years: 14.7% (male 105,763,058/female 93,903,845) 
#25-54 years: 47.2% (male 327,130,324/female 313,029,536) 
#55-64 years: 11.3% (male 77,751,100/female 75,737,968) 
#65 years and over: 9.6% (male 62,646,075/female 68,102,830) (2014 est.)


age.b <- c(0,15,25,55,65,80)

#Data for China
age.c <- rep(NA, 1000)
i <- 1
j <- i+170
age.c[i:j] <- runif(171,0,15)  
i <- j+1
j <- i+146
age.c[i:j] <- runif(147,15,25)
i <- j+1
j <- i+471
age.c[i:j] <- runif(472,25,55)
i <- j+1
j <- i+113
age.c[i:j] <- runif(114,55,65)
i <- j+1
j <- i+095
age.c[i:j] <- runif(96,65,80)
den.c <-density(age.c) # returns the density data 

#Data for Chinese internet users, from 33rd Statistical Report of Internet Development in China (CNNIC, January 2014).
#from https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwjjrvvUvInNAhVLdz4KHa6VAJsQFggdMAA&url=http%3A%2F%2Fwww1.cnnic.cn%2FIDR%2FReportDownloads%2F201404%2FU020140417607531610855.pdf&usg=AFQjCNEuKJbXON9y_nkjuvfjcjNS445cfg&sig2=gEIJ1w00TJ9tCJCFWmfVmQ
age.categories.c.i <- c(0, 10,19,29,39,49,59,80)
age.prop.c.i <- c(2.1,24.5, 30.7, 23.4, 12, 5.2, 2.1)

ld <- 1000
age.ci <- NA

for (i in 1:length(age.prop.c.i)){
  age.ci <- c(age.ci, runif(round(ld*(age.prop.c.i[i]/100)), age.categories.c.i[i], age.categories.c.i[i+1]))
}
hist(age.ci)


pdf(paste(fig.loc,"hist_hist_age.pdf", sep=""), width=13, height=8)
hist(age.c, breaks=c(0,15,25,55,65,80), col=rgb(1,0,0,.5), ylim=c(0,0.03), 
     xlab="Age", main="Histogram of Age in China (light red) and Qualtrics (green)") 
hist(age[r], breaks=c(min(age[r]),15,25,55,65,max(age[r])), col=rgb(0,1,0,.5), add=T) 
dev.off()

pdf(paste(fig.loc,"hist_den_age_comparison.pdf", sep=""), width=8, height=6)
hist(age.c, breaks=c(0,15,25,55,65,80), col=rgb(1,0,0,.5), ylim=c(0,0.05), 
     xlab="Age", main="Histogram/Density of Age in \n China, Chinese Internet, and Qualtrics") 
text(60,0.015,"Chinese", col=rgb(1,0,0,1))
#plot(den, col="cyan3", main="Age of Respondents") 
hist(age.ci, breaks=age.categories.c.i, col=rgb(0,0,1,0.5), add=TRUE)
text(10,.032, "Chinese Internet Users", col=rgb(0,0,1,1))
polygon(den, col=rgb(0,1,0,.5), border="black")
text(45,0.03,"Qualtrics", col=rgb(0,1,0,1))
rug(jitter(age, amount=0.5), col=rgb(0,1,0,0.5))
dev.off()

###Age data from the UN
un.age <- read.csv("./MASTER CODE/UNdata_age.csv")
un.age2 <- subset(un.age, Area=="Total" & Sex=="Both Sexes" & (Age %in% seq(0,99,1)))
un.age3 <- un.age2[, c("Age", "Value")]
un.age3$Age <- as.numeric(as.character(un.age3$Age))
dat <- transform(un.age3, density = Value / sum(Value))
plot(density ~ Age, data = dat, type = "n")
lines(density ~ Age, data = dat, col = "red")

pdf(paste(fig.loc,"hist_den_age_comparison2.pdf", sep=""), width=8, height=6)
plot(c(0,100),c(0,0.05), ylab="Density", xlab="Age", type="n", main="Density of Age in \n China, Chinese Internet, and Qualtrics")
text(65,0.015, "China's Age Density", col="red")
hist(age.ci, breaks=age.categories.c.i, col=rgb(0.2,1,0.2,.5), add=TRUE)
text(10,.032, "Chinese \n Internet Users", col=rgb(0.2,1,0.2))
polygon(den, col=rgb(0,0,1,0.5), border="black")
text(55,0.03, "Qualtrics' Age Density", col=rgb(0,0,1))
lines(density ~ Age, data = dat, col = "red", lwd=4)
rug(jitter(age, amount=0.5), col=rgb(0,1,0,0.5))
dev.off()

#Ethnicity
#simulating ethnicity 
d$d4 <- sample(c(1,2,99), size=nrow(d), replace=TRUE, prob=c(.42, .42, .16))
var <- d$d4
out <- rep(NA, length(var))
order.r <- rep(NA, length(var))
r <- !is.na(as.numeric(as.character(var)))
out[r] <- as.numeric(as.character(var))[r]
d$ethnicity <- out
d$ethnicity.or <- order.r
d$ethnicity2 <- as.factor(d$ethnicity)
d$ethnicity2 <- revalue(d$ethnicity2, c("1"="Han", "2"="Minority", "99"="No Answer"))

num.han <- round(sum(d$ethnicity2=="Han", na.rm=T)/sum(!is.na(d$ethnicity2)),2)
pdf(paste(fig.loc,"hist_ethnicity.pdf", sep=""), width=8, height=6)
barplot(table(d$ethnicity2)/sum(!is.na(d$ethnicity2)), main="Ethnicity")
text(0.7, 0.8, paste(num.han, "%"))
dev.off()




# D2. In what province do you live? 
#identifying. cannot be reproduced. contact the authors with questions.
region.mapping <- read.csv(file="../../../Original_data/regions.csv")

d$region <- factor(d$d2, levels=region.mapping$Code, labels=region.mapping$Region)

title <-"In what province do you live?"

reported.region <- sort(table(d$region))/length(d$region)
write.table(reported.region, file=paste(fig.loc, "reported_region.csv", sep=""))

pdf(paste(fig.loc,"region_demographics.pdf", sep=""), width=8*1.3, height=6*1.3)
par(mar=c(5,8,4,2))
barplot(reported.region, horiz=TRUE, las=1, cex.names=1, main="Reported Region")
dev.off()

### Analyzing Distribution of na1.v and na2.v

library(plot3D)

na1.v <- d$na1.v
na2.v <- d$na2.v


d$na2.v.o <-  as.factor(d$na2.v.o)
# 01 No formal education 02 Elementary school 03 Middle school
# 04 High school
# 05 College 06 Masters 07 Doctoral
d$na2.v.o <- revalue(d$na2.v.o, c("1"="Too Much", "2"="About Right", "3"="Too Little", "8"="Don't Know", "9"="Refuse to Answer"))

prop <- table(d$na2.v.o)/sum(!is.na(d$na2.v.o))
prop.other <- c(.169, .273, .36, .186, .013)
prop.data <- rbind(prop, prop.other)
row.names(prop.data) <- c("Our Sample", "RCCC 2012")
prop.data
barplot(prop.data, main="Does China rely on military strength ... (na2)", cex.names=0.7, beside=TRUE, legend=rownames(prop.data))
pdf(paste(fig.loc,"na2comp1.pdf", sep=""), width=8, height=6)
barplot(prop.data, main="Does China rely on military strength ... (na2)", cex.names=0.7, beside=TRUE, legend=rownames(prop.data))
dev.off()



prop2 <- as.data.frame(prop/sum(prop[1:3]))[1:3,]
prop.other2 <- prop.other[1:3]/sum(prop.other[1:3])
prop.data2 <- cbind(prop2, prop.other2)
colnames(prop.data2) <- c("", "Our Sample", "RCCC 2012")
prop.data2 <- t(prop.data2)
colnames(prop.data2) <- prop.data2[1,]
prop.data2 <- prop.data2[2:3,]
class(prop.data2) <- "numeric"

barplot(prop.data2, main="Does China rely on military strength ... (na2)", cex.names=0.7, beside=TRUE, legend=rownames(prop.data), args.legend=list(x="topleft"))

pdf(paste(fig.loc,"na2comp2.pdf", sep=""), width=8, height=6)
barplot(prop.data2, main="Does China rely on military strength ... (na2)", cex.names=0.7, beside=TRUE, legend=rownames(prop.data), args.legend=list(x="topleft"))
#x=5, y=.55))
dev.off()


##  Create cuts:
x_c <- cut(na1.v, 11)
y_c <- cut(na2.v, 3)

##  Calculate joint counts at cut levels:
z <- table(x_c, y_c)

##  Plot as a 3D histogram:
hist3D(z=z, border="black", xlab="\n Support for \n National Honor (na1)", ylab="\n China Should \n Use More Force (na2)", theta=-40, phi=20, main="Distribution of Pre-Scenario Responses")

p <- recordPlot()
pdf("../../Figures/na1na2.pdf", width=1.5*w, height=1.5*h)
p
dev.off()

## Analyzing aso ###

rm(list=ls())
d0 <- read.csv("./MASTER CODE/realanalysis.csv")
d1 <- read.csv("../../../Original_data/aso/150515RealHistorys-aso.csv")
varstokeep <- c("V1", "reputation",	"strength/action","sovereignty", "honor", "biding.time", "peace", "domestic", "interests", "nationalism", "deference", "complex", "nonsense.missing")
d2 <- d1[, names(d1) %in% varstokeep]
d0$V1<- d0$V1
d <- merge(d0,d2, by="V1")
remove(d0, d1, d2)

asovars <- varstokeep[-1]

asovarsindex <- names(d) %in% asovars
m <- is.na(d[, asovarsindex])
d[, asovarsindex][m] <- 0

aso.means <- colMeans(d[, asovarsindex])
aso.means <- aso.means[1:10]
aso.means.n <- as.numeric(aso.means)[1:10]

w <- 5
h <- 4

# Fitting Labels 
par(las=2) # make label text perpendicular to axis
par(mar=c(5,8,4,2)) # increase y-axis margin.
bp <- barplot(aso.means, xlab="Proportion ", horiz=TRUE, width=0.7, main="Reasons Given")
text(x=aso.means.n-0.006, y=bp, round(aso.means, 2),cex=1) 
p <- recordPlot()

pdf("../../Figures/aso1.pdf", width=1.5*w, height=1.5*h)
p
dev.off()

#-------------------------------------------------------------------------
#National Honour (Prequestions) Replication
#-------------------------------------------------------------------------

rm(list = ls())

# Customized Stargazer Function, without asterisk inflation
stargazerAD <- function(table, title="Default", filename="Figures/default.tex", dep.var.labels=NULL) {
  a <- stargazer(table, title=title,  out=filename, 
                 align=TRUE, star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c("\\dagger", "*", "**", "***"), 
                 notes="$^{\\dagger} p<0.1$; $^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$", notes.append=FALSE, single.row=TRUE)
  return(a)
}

setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

d <- read.csv("./MASTER CODE/realanalysis.csv") 

#when prequestions asked?

ggplot(d,aes(x=start.time.n,y=pre.questions))+stat_smooth(se=F) + 
  geom_vline(xintercept = d$start.time.n[1850])


df <- d[d$start.time.n<d$start.time.n[1850],]

m <- lm(asc~pre.questions, data=df)
stargazerAD(m, title="Prequestions Effect on Approval (Real History)", file="../../Figures/real_preq.tex")

#-------------------------------------------------------------------------
#Raw Results
#-------------------------------------------------------------------------
rm(list = ls())
library(ggplot2)
library(coefplot)
library(stargazer)
library(scales)
library(dplyr)
detach("package:plyr", unload=TRUE)
setwd("~/Dropbox/Dafoe-Weiss-Empirics/18-04-12-ISQ-Raluca/Code/Replication File ISQ/")

d <- read.csv("./MASTER CODE/realanalysis.csv")
full_data <- d
myvars <- names(full_data) %in% c("his", "pro", "ADIZ", "ADIZp", "eli.f", "eli.c", "asc") 
small_data <- full_data[myvars]
small_data$ID <- seq.int(nrow(small_data))

#Constructing Appropriate Comparisons 
#Treatments are assigned independently, except vague threat and defiance occurs only if vague threat
# only one of biding time or cost of war can occur 

# the control is those without biding time and those without cost of war (only one can happen)
small_data$new_eli.f[small_data$eli.f==1 ] <- 1
small_data$new_eli.f[small_data$eli.f==0 & small_data$eli.c==0] <- 0

#similarly for cost of war

small_data$new_eli.c[small_data$eli.c==1] <- 1
small_data$new_eli.c[small_data$eli.c==0 & small_data$eli.f==0] <- 0

#making a dataset for each treatment 

#history
histvars <- names(small_data) %in% c("his", "asc")
hist_data <- small_data[histvars]

hist_data$his_type[small_data$his==1] <- "Treatment"
hist_data$his_type[small_data$his==0] <- "Control"

sum(is.na(hist_data$his)) ## no missing for history treatment
sum(is.na(hist_data$asc)) ## some missing for asc. we'll take them out for the ggplot

hist_data <- na.omit(hist_data) ##taking out missing values 

hist_data2 <- hist_data %>% 
  group_by(his,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

hist_data2$overall = "History"

hist_gg <-ggplot(hist_data2, aes(x=as.factor(asc), y=perc, group=as.factor(his), fill=as.factor(his)))
hist_gg2 <- hist_gg +
  geom_bar(stat= "identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

hist_gg3 <- hist_gg2 +  
  ylab("Percentage") + 
  xlab (" ")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))




#Vague Threat

ADIZvars <- names(small_data) %in% c("ADIZ", "asc")
ADIZ_data <- small_data[ADIZvars]

ADIZ_data$ADIZ_type[small_data$ADIZ==1] <- "Treatment"
ADIZ_data$ADIZ_type[small_data$ADIZ==0] <- "Control"

sum(is.na(ADIZ_data$ADIZ_type)) ## no missing 
sum(is.na(ADIZ_data$asc)) ## some missing for asc. we'll take them out for the ggplot

ADIZ_data <- na.omit(ADIZ_data) ##taking out missing values 

ADIZ_data2 <- ADIZ_data %>% 
  group_by(ADIZ,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))


ADIZ_data2$overall = "Vague Threat"

ADIZ_gg <-ggplot(ADIZ_data2, aes(x=as.factor(asc), y=perc, group=as.factor(ADIZ), fill=as.factor(ADIZ)))
ADIZ_gg2 <- ADIZ_gg +
  geom_bar(stat="identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

ADIZ_gg3 <- ADIZ_gg2 +  
  ylab(" ") + 
  xlab (" ")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))


#Biding Time (eli.f)

elifvars <- names(small_data) %in% c("new_eli.f", "asc")
elif_data <- small_data[elifvars]

elif_data$elif_type[small_data$new_eli.f==1] <- "Treatment"
elif_data$elif_type[small_data$new_eli.f==0] <- "Control"

sum(is.na(elif_data$elif_type)) ## no missing 
sum(is.na(elif_data$asc)) ## some missing for asc. we'll take them out for the ggplot

elif_data <- na.omit(elif_data) ##taking out missing values 

elif_data2 <- elif_data %>% 
  group_by(new_eli.f,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

elif_data2$overall = "Biding Time"

elif_gg <-ggplot(elif_data2, aes(x=as.factor(asc), y=perc, group=as.factor(new_eli.f), fill=as.factor(new_eli.f)))
elif_gg2 <- elif_gg +
  geom_bar(stat="identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

elif_gg3 <- elif_gg2 +  
  ylab("Percentage") + 
  xlab ("Government Approval")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))

#Cost of war (eli.c)

elicvars <- names(small_data) %in% c("new_eli.c", "asc")
elic_data <- small_data[elicvars]

elic_data$elic_type[small_data$new_eli.c==1] <- "Treatment"
elic_data$elic_type[small_data$new_eli.c==0] <- "Control"

sum(is.na(elic_data$new_eli.c)) ## no missing 
sum(is.na(elic_data$asc)) ## some missing for asc. we'll take them out for the ggplot

elic_data <- na.omit(elic_data) ##taking out missing values 

elic_data2 <- elic_data %>% 
  group_by(new_eli.c,asc) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

elic_data2$overall = "Cost of War"

elic_gg <-ggplot(elic_data2, aes(x=as.factor(asc), y=perc, group=as.factor(new_eli.c), fill=as.factor(new_eli.c)))
elic_gg2 <- elic_gg +
  geom_bar(stat="identity", position="dodge") + 
  scale_y_continuous(labels=percent) + 
  facet_wrap(~overall)

elic_gg3 <- elic_gg2 +  
  ylab(" ") + 
  xlab ("Government Approval")  + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))   + 
  scale_x_discrete(breaks=c(1,2,3,4,5),labels= c("S Disaprove", "Disaprove", "Neither", "Approve", "S Approve")) + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name = "Treatment Status") + 
  coord_cartesian(ylim = c(0,.5))


#putting them together

library(ggpubr)

ggarrange(hist_gg3, ADIZ_gg3, elif_gg3, elic_gg3, 
          ncol = 2, nrow = 2,common.legend = TRUE, legend = "right")
ggsave(file="../../Figures/real_vertical_only4.eps")


#--------------------------------------------------------------------------------------- 
#Coefficients Plot from Regression

#Regression Plot

treat1 <- lm(asc~his, data=hist_data)
treat3 <- lm(asc~ADIZ, data=ADIZ_data)
treat5 <- lm(asc~new_eli.f, data=elif_data)
treat6 <- lm(asc~new_eli.c, data=elic_data)

mlplot <- multiplot(treat1, treat3, treat5, treat6, 
                    predictors = c("his", "ADIZ", "new_eli.f", "new_eli.c"), title="Effect on Approval, Hist", xlab="Change in Approval")

mlplot + theme_bw() + 
  scale_y_discrete(label=c("History", "Vague Threat", 
                           "Biding Time", "Cost of War")) + 
  theme(legend.position="none") + 
  scale_color_manual(values=c("black", "black", "black", "black", "black", "black"))
ggsave(file="../../Figures/real_raw_coefplot_only4.eps")

#Error Bar Plot


#History

gg1 <- ggplot(hist_data, aes(as.factor(his_type), asc, fill=as.factor(his_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("History") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  coord_cartesian(ylim=c(3.2,3.6)) + 
  theme(panel.background = element_blank(), 
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.line.x = element_blank(),
        axis.line.y = element_line(colour = "Black", linetype = "solid"))+
  scale_y_continuous(limits = c(3.2,3.6),  expand = c(0,0), oob=rescale_none) 



#Vague Threat

gg3 <- ggplot(ADIZ_data, aes(as.factor(ADIZ_type), asc, fill=as.factor(ADIZ_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("VT") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  scale_y_continuous(limits = c(3.2,3.6),  expand = c(0,0), oob=rescale_none)  + 
  coord_cartesian(ylim=c(3.2,3.6)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_blank(),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.y =element_blank(), 
        axis.text.y = element_blank(),
        axis.line.x = element_blank())


#Biding Time

gg5 <- ggplot(elif_data, aes(as.factor(elif_type), asc, fill=as.factor(elif_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("Biding") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  scale_y_continuous(limits = c(3.2,3.6),  expand = c(0,0), oob=rescale_none)  + 
  coord_cartesian(ylim=c(3.2,3.6)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_blank(),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.y =element_blank(), 
        axis.text.y = element_blank(),
        axis.line.x = element_blank())

#Cost of War 

gg6 <- ggplot(elic_data, aes(as.factor(elic_type), asc, fill=as.factor(elic_type)))+ 
  stat_summary(fun.y=mean, geom="bar", width=1, position= "dodge") +
  stat_summary(fun.data =mean_se, geom="errorbar",  width=.2, size=1, color="black",
               position=position_dodge(width = .7)) +
  ylab(" ") +   xlab("Cost of War") + 
  scale_fill_manual(labels=c("Control", "Treatment"), values=c("goldenrod1", "grey40"), name="Treatment Status") + 
  scale_y_continuous(limits = c(3.2,3.6),  expand = c(0,0), oob=rescale_none)  + 
  coord_cartesian(ylim=c(3.2,3.6)) + 
  theme(panel.background = element_blank(), 
        axis.line.y = element_blank(),
        axis.ticks.length=unit(0.3,"cm"),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), 
        axis.ticks.y =element_blank(), 
        axis.text.y = element_blank(),
        axis.line.x = element_blank())



ggarrange(gg1, gg3,  gg5, gg6, 
          ncol = 4, nrow = 1,common.legend = TRUE, legend = "right")

ggsave(file="../../Figures/real_subst_only4.eps")

#End. 